% Rounding.
%

phi = (1+sqrt(5))/2
n = 40;
f = fibonacci(1:n);
f(2:n)./f(1:n-1)
f(2:n)./f(1:n-1) - phi
format long e
n = (1:40)';
f = (phi.^(n+1) - (1-phi).^(n+1))/(2*phi-1)
f = round(f)

% check the commands fix, ceil and floor
fix(f)
ceil(f)
floor(f)
%fix(f) rounds each element of f to the nearest integer toward zero.
% For positive f, the behavior of fix is the same as floor.
% For negative f, the behavior of fix is the same as ceil.

%

% Floating point arithmetic
%floatgui

%
% Commands related with the IEEE arithmetic
%


eps %distance from 1 to the next larger floating point number = eps(1)
eps(1/2) %distance from 1/2 to the next larger floating point number
eps(2) %distance from 2 to the next larger floating point number
eps(sqrt(2)) %Meaning??? distance from sqrt(2) to the next larger floating point number
realmax %returns the largest finite floating-point number
realmin %returns the smallest positive normalized floating-point number
realmin*eps %smallest positive denormal number
0/0
inf/inf
inf - inf

%

%%%%%%%%%%%%%%%%%
realmax + 1
realmax+1-realmax
realmax-realmax+1

%

realmax == realmax +1
isequal(realmax, realmax+1)
2*realmax
realmax + 2^(1023-52)
realmax + 2^(1023-53)
realmax + 2^(1023-54)
realmin
realmin/2
realmin/2^(52)
realmin/2^(53)

%
% Single point arithmetic
%

realmax('single')
realmin('single')
log2(realmin('single'))
eps('single')
log2(eps('single'))
isequal(realmax('single'),(2-2^(-23))*2^127)

%%%%%%%%%%%%% Creating 'single-precision' variables%%%%%%%%%%
A = randn(1000); % A is double precision
As = single(A); % We ¨chop¨ A to single precision
Bs = randn(1000,'single'); % Bs is created directly in single
% Working in single precision is faster
tic;[V,D] = eig(A); toc
tic;[Vs,Ds] = eig(As); toc
% Working in single precision is less accurate
d = diag(D); ds = diag(Ds);
do = sort(d); dso=sort(ds);
max( abs(do - dso) ./ abs(do))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% Curiosities of the floating point arithmetic: Which is the right way to sum?
%

x = 1; %primero suma el 1 y los numeros pequeños (delta) los ignora, por eso da 1
delta = 1/2^(53);
for j = 1:2^(20)
    x = x + delta;
end
x

%
% Exercise: Do the same sum changing the order.
%

x = 0; %primero suma los numeros pequeños y el 1 al final, por lo que da el valor real
delta = 1/2^(53);
for j = 1:2^(20)
    x = x + delta;
end
x = 1 + x

%
% IF POSSIBLE ADD THE NUMBER WITH SMALLER ABSOLUTE VALUE FIRST
%
% Exercises:
% 1.- Operations in IEEE-Arithmetic don't fulfill the associative properties of the product and the sum. Example:

x = 1; y = 2^(-53); z = 2^(-53);
isequal(x+(y+z),(x+y)+z)

isequal(realmin*(realmin*realmax),(realmin*realmin)*realmax)
% realmin*realmax = numero que entiende y al multiplicarlo por realmin
% tambien entiende
% realmin*realmin = no entiende y al multiplicarlo por realmax, tampoco
% entiende 

% Give an example where the product associative property is not fulfilled in the IEEE arithmetic and explain why.
% 2.- The function f_1 (x,\delta) = \cos(x + \delta) - \cos(x) can be transformed into another form, f_2 (x,\delta), using the trigonometric identity:
% \cos(\phi)- cos (\varphi) =-2 \sin (\frac{\phi +\varphi}{2} ) \sin(\frac{\phi -\varphi}{2})
% Thus, f_1 and f_2 have the same values, in exact arithmetic, for any given argument values x and \delta.
% Write a MATLAB script that will calculate g_1(x,\delta)=f_1(x,\delta)/\delta+\sin(x) and $g_2(x,\delta)= f_2(x,\delta)/\delta+\sin(x) for x = 3 and \delta = 1.0e-11. Explain the difference in the results of the two calculations.

function [f1] = func1(x,delta)
f1 = cos(x + delta) - cos(x);
end

function [f2] = func2(x,delta)
f2 = (-2*sin((2*x + delta)/2))*(sin(delta/2));
end


% 3.- The function f_1(x_0,h)= \sin(x_0 + h)-\sin(x_0) can be transformed to another form, f_2 (x_0,h), using the trigonometric formula:
% \sin(\phi)- \sin(\varphi) =2 \cos(\frac{\phi +\varphi}{2} ) \sin(\frac{\phi -\varphi}{2})
% Thus, f_1 and f_2 have the same values, in exact arithmetic, for any pair of arguments x_0 and h.
% Find a formula that avoids cancellation errors for computing the approximation \frac{f(x_0 + h)-f(x_0)}{h} to the derivative of f(x)= \sin(x) at x = x_0 . Write a MATLAB  program that implements this formula and computes an approximation of f'(.5), for h=1e-20,\cdots, 1. Compare your results to those in Example 3.1. Explain the difference between your results and the results reported in Example 3.1.
% 4.- The number \frac{3}{8} = 2.666... has no exact representation in any decimal floating-point system with finite precision t. Is there a finite floating-point system in which this number does have an exact representation? If yes, then describe one such system. Answer this question for the irrational number \pi
% 5.- The roots of the quadratic equation x^2 -2bx +c=0 with b^2 >c are given by:
% x_{1, 2} = b \pm \sqrt{b^2-c}
% Remember that x_1x_2=c. The MATLAB scripts calculate these roots in two different ways:
% (a)      \matrix{ x_1 = b + sqrt(b ^2-c)  \cr x_2 = b - sqrt(b^2-c) }
% (b)   if b > 0
%             x_1 = b + sqrt(b ^2-c) ;
%             x_2 =c / x_1;
%         else
%             x_2 = b - sqrt(b^2-c) ;
%             x_1 =c/x_2;
%         end
% Which algorithm gives a more accurate result in general? Try to answer this question without any computing.
% 6.- Write a MATLAB program that will:
%     (a)   Sum up \frac{1}{n}  for n = 1,2,...,100,000
%     (b)  Round each number - to four decimal digits, and then sum them up in five-digit decimal arithmetic for      n = 1,2,...,100,000.
%     (c) Sum up the same rounded numbers (in four-digit decimal arithmetic) in reverse order, i.e., for n = 100,000,...,2, 1. Compare the three results and explain your observations.
% 7.- In statistics, one often needs to compute the quantities \bar{x}=\frac{1}{n}\sum_{i=1}^n x_i and s^2 = \frac{1}{n} \sum_{i=1}^n (x_i-\bar{x})^2, where x_1,x_2,...,x_n,  are the given data. It is easy to see that s^2 can also be written as s^2=\frac{1}{n}\sum_{i=1}^{n} x_i^2-\bar{x}^2.
%     (a) Which of the two methods to calculate s^2 requires fewer operations? Assume that \bar{x} has already been calculated when counting the operations.
%     (b) Which of the two methods is expected to give more accurate results for s^2 in general?
%     (c) Write a MATLAB script to check whether your answer to the previous question was correct.
%
% This one is more theoretical (use pena and paper).
% 8.- With the exact rounding of Section 3.5, each elementary operation has a relative error that is bounded in terms of the rounding unit, which is, in MATLAB, eps. Now consider exponentiation, which is implemented as x^y = e^{y\ln(x)}. Estimate the relative error in calculating x^y in floating point, assuming fl(x + y) = (x + y)(1 + \epsilon), | \epsilon | \leq eps, and that everything else is exact. Show that the bounds for elementary operations do not hold for exponentiation when x is very large.